####################################################################################
## xu xian feng
## 2022/03/30
## STAD analysis
####################################################################################

source ~/20220915_gastric_multiple/dna_combinePublic/config/config.sh
export config_file=~/20220915_gastric_multiple/dna_combinePublic/config/config.sh

####################################################################################
## 第一部分
## NMU样本的质控及突变CCF的计算
####################################################################################
##########################################
#### 得到NMU的MSS数据信息
sh ${scripts_path}/module/pipeline_nmu.sh
## Mutation的按照是否与胃癌共享定义为超早期
${Rscript} ${scripts_path}/plot/MutationTime_NewTime.R \
--ccf_file ${MutationTime_path}/result/All_CCF_mutTime.tsv \
--sample_info ${config_path}/tumor_normal.class.list \
--out_path ${MutationTime_path}/result

## NMU的MSI数据整理
sh ${scripts_path}/module/pipeline_nmu_msi.sh

## 所有最后分析的所有样本列表
## 包含了MSI
cat ${config_path}/tumor_normal.class.list > ${config_path}/tumor_normal.class.MSS_MSI.list
cat ${config_path}/tumor_normal.class.MSI.list | grep -v  -w "ID" >> ${config_path}/tumor_normal.class.MSS_MSI.list

##########################################
#### 整理NMU基线
## 包含MSI的样本
## 整理年龄、性别、吸烟、饮酒等信息
## 注释整体和CDS的覆盖区域、倍体、纯度、质控信息

## 所有NMU样本的基线信息（STAD-useCombine.Patient.tsv）
## 按照IM+IGC,IM+DGC,IM+IGC+DGC，分成三类比较（STAD-useCombine.BaseLineCompare.csv）
${Rscript_mutationTime} ${scripts_path}/CompareBaseLine.R \
--tumor_list ${config_path}/tumor_normal.class.MSS_MSI.list \
--baseline_file ${config_path}/STAD_MutipleReigon_baseline.addAlcoholFreq.tsv \
--qc_file ${Qc_path}/Summary_Qc.tsv \
--burden_all_file ${Qc_path}/Burden.coverage10x.Autosomal.txt \
--burden_cds_file ${Qc_path}/Burden.coverage10x.Autosomal.cds.txt \
--purity_file ${Titan_path}/Purity_titan.final.tsv \
--out_dir ${config_path}

## 在每个样本的基线基础上，增加突变负荷信息（mutBurden.tsv）
${Rscript} ${scripts_path}/mutBurden/mutBurden_plot.addMSI.R \
--info_file ${config_path}/STAD-useCombine.Sample.tsv \
--maf_file ${maf_path}/All_ForMutBurden.extract.maf \
--maf_msi_file ${maf_path}/All_ForMutBurden.extract.MSI.maf \
--images_path ${Images_path}/mutBurden

## 软链接
ln -snf ${Images_path}/mutBurden/mutBurden.tsv ${baseTable_path}/STAD_Info.addBurden.MSI_MSS.tsv

## 非MSI的患者其突变负荷情况
msi_high_sample=`cat ${Qc_path}/Summary_Qc.tsv | grep -v Tumor | awk -F'\t' '{if($8 > 10)print $2}' | sort -u |\
tr '\n' '|' | sed 's/|$//'`
msi_high_sample=`echo "${msi_high_sample}|JRGC00009"`
cat ${baseTable_path}/STAD_Info.addBurden.MSI_MSS.tsv | grep -v -E -w ${msi_high_sample} \
> ${baseTable_path}/STAD_Info.addBurden.MSS.tsv

####################################################################################
## GS亚型定义
## 用于分子分型
## NMU所有样本，包含了MSI的样本
sh ${scripts_path}/gistic2/gistic2.CIN.sh
## 判断CNV的亚型
${Rscript} ${scripts_path}/gistic2/classify_CIN.R \
--sample_file ${config_path}/tumor_normal.class.list \
--burden_file ${work_dir}/baseTable/STAD_Info.addBurden.MSI_MSS.tsv \
--gistic_file ${work_dir}/gistic/all_lesions.conf_99.txt \
--out_file ${work_dir}/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv



####################################################################################
## 第二部分
## 公共数据样本的突变数据整理
####################################################################################
##########################################
#### 整理公共数据库的基线
## 以下文件均在${work_dir}/public_ref/combine文件夹下面
## 整合公共数据，所有GC样本的基线和TMB（MutationInfo.combine.tsv）
${Rscript} ${scripts_path}/comparePublicData/CombinePublicData.R \
--njmu_file ${work_dir}/baseTable/STAD_Info.addBurden.MSI_MSS.tsv \
--tcga_file ${work_dir}/public_ref/TCGA/TCGA_STAD.TMB.tsv \
--tcga_clinic_file ${work_dir}/public_ref/TCGA/clinical_PANCAN_patient_with_followup.tsv \
--tcga_msi_file  ${work_dir}/public_ref/TCGA/TCGA_STAD_msi.TMB.tsv \
--oncoSG_file ${work_dir}/public_ref/OncoSG/OncoSG_STAD.TMB.tsv \
--oncoSG_clinic_file ${work_dir}/public_ref/OncoSG/OncoSG_STAD.followup.tsv \
--oncoSG_msi_file ${work_dir}/public_ref/OncoSG/OncoSG_STAD_msi.TMB.tsv \
--tmucih_file ${work_dir}/public_ref/TMUCIH/TMUCIH_STAD.TMB.tsv \
--tmucih_msi_file ${work_dir}/public_ref/TMUCIH/TMUCIH_STAD_msi.TMB.tsv \
--utokyo_file ${work_dir}/public_ref/utokyo/utokyo_STAD.TMB.tsv \
--out_path ${work_dir}/public_ref/combine

## 增加CIN的分型（MutationInfo.combine.addMolecularSubType.tsv）
${Rscript} ${scripts_path}/comparePublicData/CombinePublicData.addCNVType.R \
--input_file ${work_dir}/public_ref/combine/MutationInfo.combine.tsv \
--njmu_clinic_file ${work_dir}/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv \
--tcga_clinic_file ${work_dir}/public_ref/TCGA/stad_tcga_pan_can_atlas_2018_clinical_data.tsv \
--oncoSG_clinic_file ${work_dir}/public_ref/OncoSG/stad_oncosg_2018_clinical_data.tsv \
--tmucih_clinic_file ${work_dir}/public_ref/TMUCIH/egc_tmucih_2015_clinical_data.tsv \
--out_path ${work_dir}/public_ref/combine

## 标记样本人种（MutationInfo.combine.addMolecularSubType.Race.tsv）
${Rscript} ${scripts_path}/plot/mutBurden.population.R \
--input_file ${work_dir}/public_ref/combine/MutationInfo.combine.addMolecularSubType.tsv \
--tcga_file ${work_dir}/public_ref/TCGA/stad_tcga_pan_can_atlas_2018_clinical_data.tsv \
--out_path ${work_dir}/public_ref/combine

##################################### 
#### 合并公共数据库的maf文件
## 得到最后使用的maf列表
## 同一人同一类型样本突变合并
## 用于突变率的估计和显著突变基因的计算
## 去除IGC+DGC的样本（NMU的5个样本），只考虑一种肿瘤类型的
## All_use.IM.maf
## All_use.maf
sh ${scripts_path}/module/getmaf_public_unique.sh

## MSI的
sh ${scripts_path}/module/getmaf_public_unique.msi.sh

## 同一人突变保留
## All_use.addVAF.maf
sh ${scripts_path}/module/getmaf_public_allSample.sh
## All_use.addVAF.maf和All_use.maf的主要差别在于
## 1、All_use.addVAF.maf增加了t_alt_count和t_ref_count
## 2、All_use.addVAF.maf中NMU的样本为原始的ID号，同一个人多个样本保留，IM的样本也在里面，IGC+DGC的样本也在里面；
## 3、All_use.maf中为NMU的GC样本，无IGC+DGC的5个人，All_use.IM.maf中为所有人IM样本
## MSI的患者
## All_use.addVAF.MSI.MSI.maf，为MSI的患者
sh ${scripts_path}/module/getmaf_public_allSample.msi.sh


##################################### 
#### CNV的除了NMU的只有TCGA存在，整理
sh ${scripts_path}/module/getCNV_public.sh


##################################### 
#### 合并TCGA的RNA
## 1、提取TCGA明确laruen分型的样本
## 2、与NJMU所有样本去除批次效应
## CombineTMM.DNAUse.NJMU_TCGA.MergeMutiSample.tsv，合并NJMU的多样本，用于分析表达变化
## CombineTMM.DNAUse.NJMU_TCGA.tsv，提取用到的NJMU样本，用于分析突变情况
## TCGA.FilterLowExpression.TMM.tsv，TCGA单独的表达矩阵
## NJMU.FilterLowExpression.TMM.tsv，NMU单独的表达矩阵
## NJMU.FilterLowExpression.MergeMutiSample.TMM.tsv，NMU单独的表达矩阵
sh ${scripts_path}/module/getRNAData.sh


####################################################################################
## 第三部分，驱动基因相关鉴定和突变率计算
####################################################################################
##################################### 
#### 所有基因突变率的计算，IM、IGC、DGC、GC
## MutRate.tsv
${Rscript} ${scripts_path}/plot/mutRate_compute.R \
--maf_cancer_file ${maf_public_path}/All_use.maf \
--maf_im_file ${maf_public_path}/All_use.IM.maf \
--images_path ${Images_path}/mutRate \
--info_file ${work_dir}/public_ref/combine/MutationInfo.combine.tsv


##################################### 
## 驱动基因的鉴定，Mutsig2CV
## 分IGC、DGC、IM分别鉴定显著突变基因
for class in ${class_type[@]}
do
echo $class
sh ${scripts_path}/mutsig/Mutsig_3_mutsigcv.sh ${class}
done


## 检查显著驱动基因的突变率在不同来源样本的突变情况
## ${mutsig_check_path}/driver.summary.tsv
${Rscript} ${scripts_path}/mutsig/Mutsig_combineInfo.R \
--sig_igc_file ${MutsigOut_path}/IGC/sig_genes.txt \
--sig_dgc_file ${MutsigOut_path}/DGC/sig_genes.txt \
--sig_im_file ${MutsigOut_path}/IM/sig_genes.txt \
--reprot_file ${work_dir}/public_ref/SMG_sort.list \
--mutRate_file ${Images_path}/mutRate/MutRate.tsv \
--out_path ${mutsig_check_path}


## 人工检查去除不可靠的基因
## NMU的检查bam文件
# ${mutsig_check_path}

## 可靠的基因列表
# ${mutsig_check_path}/cancer/cancer.smg.list
# ${mutsig_check_path}/precancer/precancer.smg.list
## 合成一个文件
echo "Gene_Symbol" > ${mutsig_check_path}/smg.list
cat ${mutsig_check_path}/igc_smg.list ${mutsig_check_path}/dgc_smg.list ${mutsig_check_path}/im_smg.list | \
grep -v Gene_Symbol | sort -u \
>> ${mutsig_check_path}/smg.list


####################################################################################
## 第四部分，构建系统发生树
####################################################################################
#### 计算每个人样本的数量
## 2个人构建进化树纳入所有突变不考虑vaf
## 超过3个人构建进化树需考虑vaf使树更可靠
${Rscript} ${scripts_path}/tree/computSampleNum.R \
--info_file ${config_path}/tumor_normal.class.MSS_MSI.list \
--out_file ${tree_path}/tumor_normal.class.forTree.list

## 描述每个样本的driver基因VAF情况以及对应样本的突变
## 为后期检查进化树
${Rscript} ${scripts_path}/tree/GetDriverEverySample.R \
--ccf_file ${MutationTime_path}/result/All_CCF_mutTime.tsv \
--ccf_msi_file ${MutationTime_path}/result/All_CCF_mutTime.MSI.tsv \
--gene_list ${work_dir}/public_ref/importTantGene.list \
--sample_info ${config_path}/tumor_normal.class.MSS_MSI.list \
--out_path ${Images_path}/DriverEverySample \
--class_order_file ${config_path}/Class_order_sub.list 

##################################### 
#### 环境配置
## 系统发生树
## 环境在950上
## 拷过去操作
export PROOT_NO_SECCOMP=1
HOME=/public/user/xxf2019

${HOME}/udocker run -v /public:/public treeomics bash

##################################### 
#### 以下操作在udocker里面进行
Xvfb -ac :11 -screen 0 1280x1024x8 &
export  DISPLAY=:11
export HOME=/public/user/xxf2019
export work_dir=~/20220915_gastric_multiple/dna_combine
source ${work_dir}/config/config.sh

##################################### 
#### 连接需要展示的基因集
## 驱动基因标注后期AI重新标记
mv /treeomics/input/SMG_sort.txt /treeomics/input/SMG_sort.bk.txt
cp -rf ${work_dir}/images/selectGCClone/GCClone_gene.all_record.list /treeomics/input/SMG_sort.txt
echo "BMP6" >> /treeomics/input/SMG_sort.txt
echo "MUC6" >> /treeomics/input/SMG_sort.txt
echo "CFTR" >> /treeomics/input/SMG_sort.txt

##################################### 跑进化树
HOME=/
## conda install -c bioconda pyensembl
## pyensembl install --release 75 --species homo_sapiens
## 使用纳入所有突变的配置文件
## 当样本数量小于3个时
## 构建进化树，纳入所有突变
cp -f /root/treeomics/treeomics/patient.py.revise /root/treeomics/treeomics/patient.py

export min_vaf=0.0001
export error=0.00000001
cat ${tree_path}/tumor_normal.class.forTree.list | grep -v Normal | awk -F'\t' '{if($8<=2)print $2}' | sort -u | xargs -P 10 -I Normal sh -c '
echo Normal
sh ${scripts_path}/tree/Treeomics_Tree.sh Normal ${config_path} ${min_vaf} ${error}
'

## 当含有大于2个样本的时候，纳入所有突变会造成进化树构建不准确，提高VAF判定为假的阈值
cp -f /root/treeomics/treeomics/patient.py.raw /root/treeomics/treeomics/patient.py 
export min_vaf=0.0001
export error=0.00000001
cat ${tree_path}/tumor_normal.class.forTree.list | grep -v Normal | awk -F'\t' '{if($8>2)print $2}' | sort -u  | xargs -P 10 -I Normal sh -c '
echo Normal
sh ${scripts_path}/tree/Treeomics_Tree.sh Normal ${config_path} ${min_vaf} ${error}
'

##################################### 拷贝文件
rm -rf ${tree_path}/Tree_file
mkdir -p ${tree_path}/Tree_file

cat ${tree_path}/tumor_normal.class.forTree.list | grep -v Normal | awk -F'\t' '{print $2}' | \
sort -u | xargs -P 10 -I Normal sh -c '
echo Normal
sh ${scripts_path}/tree/Treeomics_moveTree.sh Normal ${config_path}
'

##################################### 
#### 进化树完成人工检查以后，重新链接
for id in `cat ${Images_path}/ITH/DriverClass.tsv | grep -v Normal | awk '{print $1}'`
do
echo $id
normal=`cat ${config_path}/tumor_normal.class.list | grep ${id} | awk '{print $2}' | sort -u`
tree_class=`cat ${Images_path}/ITH/DriverClass.tsv | grep -w ${id} | awk -F'\t' '{print $3}'`
## 驱动基因的CCF变化
ln -snf ${Images_path}/DriverEverySample/${id}_Driver.pdf ${tree_path}/TreeClass_Revise/${tree_class}/${id}_Driver.pdf
done


####################################################################################
## conda导出用到的运行环境
####################################################################################

mkdir -p ${work_dir}/conda_env
cd ${work_dir}/conda_env

#### 画图和大部分统计分析用的环境
## R 4.1.0
conda env export -n test_new > test_new.yml

#### 单细胞用的scissor的环境
## R 4.1.3
## Scissor_2.0.0
## Seurat_4.2.0
conda env export -n singleCell > singleCell.yml

#### 计算CCF
## R 4.1.0
## CHAT_1.0
conda env export -n xxf_chat_new > xxf_chat_new.yml

#### 计算突变克隆时期以及产生基线表格
## R 4.1.2
## MutationTimeR_1.00.2
## tableone_0.13.2
conda env export -n mutationTimeR > mutationTimeR.yml

#### 计算受体-配体互作
## Python 3.8.16
## ktplotspy==0.1.10
## R 4.2.2
## Seurat_4.3.0
conda env export -n cpdb4.0 > cpdb4.0.yml

#### 计算进展时间
## R 4.2.0
## ggmcmc_1.5.1.1
## rjags_4-14
conda env export -n evolutionTime > evolutionTime.yml


####################################################################################
## 备份到GSA的测序数据(该ID不满足需要，重新编码)
sh ${scripts_path}/Datalist/pipeline_list_dna.sh
sh ${scripts_path}/Datalist/pipeline_list_rna.sh


cd ${work_dir}/public_data/
ln -snf NJGC_Somatic.datalist.csv NJGC_Somatic.datalist.dna.csv

## 基线列表,替换ID
cat ${config_path}/STAD-useCombine.Patient.tsv > ${work_dir}/public_data/NJGC_Somatic.baseinfo.tsv
for line in `cat ${config_path}/plotID.list | awk -F'\t' '{OFS=","}{print $1,$5}' `
do
rawid=`echo ${line} | awk -F, '{print $1}' `
plotid=`echo ${line} | awk -F, '{print $2}' `
sed -i "s/${rawid}/${plotid}/" ${work_dir}/public_data/NJGC_Somatic.baseinfo.tsv
done

## 生存fastq文件的md5
cd ${work_dir}/public_data/wes/fastq
md5sum * > ${work_dir}/public_data/wes/fastq.md5 &

cd ${work_dir}/public_data/rna/fastq
md5sum * > ${work_dir}/public_data/rna/fastq.md5 &

####################################################################################
## 备份到GSA的测序数据,重新编码
sh ${scripts_path}/Datalist/pipeline_list_dna.v2.sh
sh ${scripts_path}/Datalist/pipeline_list_rna.v2.sh

cd ${work_dir}/public_data/
ln -snf NJGC_Somatic.datalist.csv NJGC_Somatic.datalist.dna.csv

## 基线列表,替换ID
cat ${config_path}/STAD-useCombine.Patient.tsv > ${work_dir}/public_data/NJGC_Somatic.baseinfo.tsv
for line in `cat ${config_path}/plotID.list | awk -F'\t' '{OFS=","}{print $1,$4}' `
do
rawid=`echo ${line} | awk -F, '{print $1}' `
plotid=`echo ${line} | awk -F, '{print "NJ"$2}' `
sed -i "s/${rawid}/${plotid}/" ${work_dir}/public_data/NJGC_Somatic.baseinfo.tsv
done

## 生存fastq文件的md5
cd ${work_dir}/public_data/wes/fastq
md5sum * > ${work_dir}/public_data/wes/fastq.md5 &

cd ${work_dir}/public_data/rna/fastq
md5sum * > ${work_dir}/public_data/rna/fastq.md5 &



####################################################################################
## 20241019补充突变负荷(基于766个基于列表，来源于2018年的cancer cell)
## 
export gene_list_file=~/20220915_gastric_multiple/dna_combinePublic/public_ref/2018_CancerCell.targetgene.list
## 提取基因所在区域
cat ~/ref/PCAWG_Elements/web_hg19/gc19_pc.cds.use.nochr.bed | grep -w -f ${gene_list_file} | grep -v "-" |\
awk -F'\t' '{OFS="\t"}{print $1,$2,$3}' \
> ~/20220915_gastric_multiple/dna_combinePublic/public_ref/2018_CancerCell.targetgene.bed

export region_type=target766
export capture_bed=~/20220915_gastric_multiple/dna_combinePublic/public_ref/2018_CancerCell.targetgene.bed

## 靶向区域的覆盖深度
for sample in `cat ${config_path}/tumor_normal.class.MSS_MSI.list | awk '{print $3}'` 
do
echo ${sample}
echo " sh ~/dna_preprocess_scripts/mosdepth_run_region.sh ${sample} ${config_path} ${capture_bed} ${region_type}" | \
qsub -N ${sample}_"Depth" -l nodes=1:ppn=2,mem=10gb,walltime=240:00:00 -q batch -d ${Qsub_log_path}
done

## 合并深度到一个文件
echo -e "Tumor\tcoverage" > ${qc_path}/Burden.coverage10x.Autosomal.${region_type}.txt
for sample in `cat ${config_path}/tumor_normal.class.MSS_MSI.list | awk '{print $3}'` 
do
echo $sample
region=`zcat ${qc_path}/mosdepth_${region_type}/${sample}.thresholds.bed.gz | awk -F'\t' '{if($1!="X" && $1!="Y")print}' | awk -F'\t' '{sum+=$6}END{print sum}'`
echo -e ${sample}"\t"${region} >> ${qc_path}/Burden.coverage10x.Autosomal.${region_type}.txt
done

## 计算靶向基因突变的数量
echo -e "sample\tmutNum" > ${qc_path}/${region_type}.mutNum.txt
export Variant_Types="Missense_Mutation|Nonsense_Mutation|Frame_Shift_Ins|Frame_Shift_Del|In_Frame_Ins|In_Frame_Del|Splice_Site|Nonstop_Mutation"
for sample in `cat ${config_path}/tumor_normal.class.MSS_MSI.list | grep -v Tumor | awk '{print $3}'` 
do
echo ${sample}
mut_num=`cat ${maf_path}/${sample}_*_QC.maf | awk -F'\t' '{OFS="\t"}{print $1,$9}' | grep -E -w ${Variant_Types} | \
grep -w -f ${gene_list_file} | awk -F'\t' '{OFS="\t"}{if($1!~"-"){print}}' | wc -l`
echo -e "${sample}\t${mut_num}" >> ${qc_path}/${region_type}.mutNum.txt
done

## 计算突变负荷
${Rscript} ${scripts_path}/mutBurden/mutBurden_plot.addMSI.addTarget766.R \
--info_file ${work_dir}/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv \
--coverage_file ${qc_path}/Burden.coverage10x.Autosomal.${region_type}.txt \
--mutnum_file ${qc_path}/${region_type}.mutNum.txt \
--images_path ${work_dir}/baseTable/addTarget766Burden 
